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SUMMARY 

A 

A  method  of  finding  the  flexural  and  torsional  normal  modes  of  structures 
which  have  straight  stiffness  axes  and  which  are  mounted  as  cantilevers  is  given. 
The  Lanczos  method  of  minimised  iterations  is  used  to  obtain  intermediate  modes 
using  the  integro-dif ferential  equations  and  the  mass  and  stiffness  distributions. 
The  inertia  matrix  appropriate  to  the  intermediate  modes  is  tridiagonal  and  the 
stiffness  matrix  is  unit.  The  dominant  eigenvalues  and  vectors  of  the  inertia 
matrix  give  good  approximations  to  the  graver  normal  modes.  Results  of  test 
calculations  using  a  computer  program  which  also  allows  for  the  presence  of 
discrete  masses  are  given.  ^ —  > 
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1  INTRODUCTION 

In  the  early  stages  of  the  design  of  wings  of  moderate-to-large  aspect 
ratios  their  structural  properties  are  often  given  in  terms  of  spanwise  mass  and 
stiffness  distributions,  the  stiffness  properties  being  taken  to  be  concentrated 
at  a  straight  spanwise  stiffness  axis.  The  normal  modes  of  such  wings,  mounted 
as  cantilevers,  are  needed  for  investigations  of  their  oscillatory  characteris¬ 
tics  and  aeroelastic  behaviour  and  a  method  is  described  herein  whereby  these 
modes  can  be  obtained  directly  from  the  distributions  of  mass  and  stiffness  with¬ 
out  the  intermediate  step  of  either  replacing  the  continuous  mass  of  the  beam  by 
discrete  masses  and  moments  of  inertia  and  the  stiffness  by  a  flexibility  matrix 
with  a  consequent  loss  of  accuracy  or  choosing  an  adequate  set  of  assumed 
modes . 

The  method  is  based  on  the  numerical  solution  of  the  integro-dif ferential 
equations  describing  flexural  and  torsional  oscillations  (equation  (6)).  With¬ 
out  the  help  of  a  digital  computer  such  solutions  are  tedious,  calling  as  they 
do  for  the  evaluation  of  numerous  integrals.  A  digital  computer  removes  the 
tedium;  and  accuracy  can  be  achieved  by  integrating  for  a  large  number  of  span- 
wise  stations  even  if  the  method  of  numerical  integration  used  is  the  simplest. 
The  direct  iterative  solution  of  the  equations,  Stodola's  method*,  breaks  down 
when  the  frequencies  of  two  or  more  of  the  modes  are  close  together  in  just  the 
same  way  as  does  its  equivalent,  the  power  method  of  finding  the  eigenvalues  of 
matrices.  A  possible  way  of  avoiding  this,  whilst  retaining  the  valuable 
properties  of  Stodola's  method,  is  to  find  a  set  of  intermediate,  modes  which  are 
orthogonal  to  each  other  with  respect  to  stiffness  by  Stodola  iteration  and  then 
to  solve  the  eigenvalue  equation  which  has  the  inertia  and  stiffness  matrices  of 
these  modes  as  coefficients.  Not  onlj  will  the  stiffness  matrix  be  diagonal  but 

the  inertia  matrix  will  be  tridiagonal  for  this  is  the  Lanczos  method  of 
...  2 

minimised  iterations  .  For  convenience  the  intermediate  modes  will  be  referred 
to  as  Lanczos  modes.  They  can  be  looked  upon  as  a  set  of  assumed  modes  with 
special  properties. 

The  first  Lanczos  modes  are  composed  predominantly  of  the  graver  normal 
modes  just  as  the  first  mode  found  by  Stodola's  method  is  the  fundamental. 

In  fact  Lanczos  noticed  that,  if  the  frequencies  of  the  normal  modes  are  well- 
separated,  the  inclusion  of  an  extra  Lanczos  mode  adds  one  more  normal  mode  to 
the  solution  and  corrects  the  graver  normal  modes  by  only  small  amounts.  In 
general  the  frequencies  are  not  as  well-separated  as  this  but  the  economy  that 
can  be  exercised  in  the  number  of  Lanczos  modes  used  is  the  principal  attraction 
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of  the  method.  Were  it  not  possible  to  get  good  results  with  a  restricted  number 

of  Lanczos  modes  the  method  would  not  be  used  because  as  the  number  of  Lanczos 

modes  calculated  is  increased  cancellation  results  in  a  crippling  loaa  of 
3 

accuracy  . 

In  the  following  sections  an  iterative  procedure  for  transforming  the 

integro-dif ferential  equations  into  a  tridiagonal -matrix  equation  is  derived 

from  an  established  procedure  for  transforming  the  full-matrix  eigenvalue 

equation  and  some  details  of  the  computer  program  that  has  been  written  are  given. 

Wings  often  carry  discrete  masses  and,  because  these  cannot  be  represented 

accurately  in  the  simple  methods  of  integration  used,  provision  is  made  in  the 

program  for  the  inclusion  of  discrete-load  modes.  These  are  particularly  useful 

.  4 

in  investigations  m  which  the  sizes  of  the  discrete  masses  are  varied  .  Finally 
the  results  of  tests  to  establish  the  range  of  applicability  of  the  program  are 
given. 

2  THE  ITERATIVE  PROCEDURE 

2. 1  The  equivalent  matrix  equation 

The  basis  of  the  method,  described  in  section  2.2  below,  is  derived  from 
that  used  for  finding  the  dominant  eigenvalues  of  the  equivalent  matrix  equation 
which  can  be  written 

(A  -  AE)qa  =  0  (1) 


where  A  and  E  are  square  symmetric  matrices  of  structural  inertia  and  stiff- 

-2 

ness  coefficients,  A  =  w  ,  u  is  the  frequency  of  a  normal  mode  and  q^  is 
its  modal  vector. 


The  Lanczos  method  of  minimised  iteration  can  be  applied  to  find  n  inter¬ 
mediate  modes  by  the  following  sequence  of  operations'*: 

A 

(1)  Take  a  random  vector,  qj  ,  and  put  i  =  1 

(2)  Put  8?  .  =  q?Eq.  ; 

i-l  ni  _j 

(3)  Normalise  by  putting  q^  ■  8._|q.  J 

(4)  Find  the  direct  inertia  coefficient,  a.  = 

(5)  Iterate  to  obtain  next  raw  vector,  q^+1  ■ 

(6)  Orthogonalise  with  respect  to  stiffness, 
where  gQq0  =  0  ; 

(7)  Add  1  to  i  and  repeat  from  (2)  until  required  number, 
have  been  obtained  when  stop  at  (4) . 


; 


E",Aqi  ; 

,  Of  q. 


>  (2) 


A  dash  in  the  above  signifies  transposition 
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It  has  been  shown 


2,3,5 


that 


Q'AQ  =  A  (3) 

where  Q  =  (q, ,  q2>  . . . ,  qR)  , 

and  A  is  a  symmetric  tridiagonal  matrix  with  leading  diagonal 

(cij.o^,  . . . ,  a  )  and  sub-  and  super-diagonals  each  (3  j »  ®  2 »  ****®n-)^‘ 

In  the  sequence  (2)  the  have  been  normalised  so  that  Q'EQ  =  I  and 

so  the  solutions  of 

(A  -  AI)qb  =  0  (4) 

are  approximations  to  the  dominant  eigenvalues  of  equation  (1).  The  eigen¬ 
vectors,  q  ,  are  then  given  by 

qa  =  Qqb  .  (5) 

The  order,  n  ,  of  the  matrices  in  equation  (4)  need  only  be  about  twice  the 
number  of  eigenvalues  required. 

2.2  The  integro-dif ferential  equations 

The  integro-dif ferential  equation  which  governs  a  wing  vibrating  in 
flexure  and  torsion  can  be  written 


e"(n) 

2 

=  0) 

*«(n)*ce(n) 

"1  l 

j  dnj  1  |u(n2>C(n2)  +  u? (n2) © (n2> f dri2 
n  n  | 

e'(n) 

*0!;(n>*ee(n) 

} 

1  {y£;(n2)c(n2)  +  u*2  (n2>  ©  (n2>  f  dn2 

am  _ 

n  _ 

where  n  is  a  spanwise  coordinate,  zero  at  the  root  and  unity  at  the  tip, 

4(n)  is  a  scaled  heave  displacement, 

©Cm)  is  an  incidence  displacement, 

<t>ge(n),  4>^g(h),  are  the  direct  and  cross-flexibilities, 

y(n)  is  the  line  density  of  the  mass, 

p£(n)  is  the  line  density  of  the  first  moment  of  mass  about  the  stiffness 

axis, 

__2 

ytc  (n)  is  the  line  density  of  the  second  moment  of  mass  about  the  local 
centre  of  gravity, 

and  dashes  here  signify  differentiation  with  respect  to  n  . 
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The  2><2  matrix  of  flexibilities  includes  crossflexibilities  so  that 
account  may  be  taken  of  the  special  properties  of  fibre-reinforced  plastics  with 
asymmetric  lay-ups. 

The  essentials  of  the  method  can  be  demonstrated  by  considering  only  heave 
and  so  the  equation  considered  is 


=  a) 

where  is  the  bending  moment  distribution,  the  repeated  integral  in 

equation  (6).  For  our  immediate  purposes  the  suffices  t,  can  be  omitted. 

In  the  Stodola  method  the  sequence  of  operations  is: 

(1)  Choose  a  deflection  mode  £  containing  components  of  all  the  graver 
normal  modes  and  put  i  =  1  ; 


(7) 


(2) 


Put  =  C  ; 


(3)  Make  the  mode  orthogonal  to  previous  modes  with  respect  to  mass  by 

...  ...  i-1 

putting  cH  =  +  23  ^k^k  w^ere  t*le  <jc  are  chosen  to  make 


k=l 


J  =  0  (i  >  1  only); 

rv  1  K- 


1 


l 


(4) 

(5) 

(6) 


'1 


Find  the  bending  moment  distribution,  *  J  dn^  /  dr^  i 

n 

Find  the  curvature  distribution,  +  =  4> ij;  f ^  ; 

.,(j  +  0 


>  (8) 


Integrate 


twice  to  obtain  the  deflection  mode  Z 


(j  +  1) 


and  normalise  by  putting  +  +  ,  where 

is  determined  by  the  condition  +  =  1  » 

1  'tip 

j  2 

(7)  Test  for  convergence.  If  /  +  dp  >  e  ,  a  suitable 

0V1  1 

small  number  related  to  the  numerical  accuracy  of  the  computer, 
add  1  to  j  and  go  to  (3);  otherwise,  accept  the  normal  mode, 

+  ,  and  natural  frequency,  w.  ■  +  ,  add  1  to  i 

and  go  to  (2). 
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T 


(!) 

(2) 

(3) 

(4) 

(5) 

(6) 


In  the  proposed  method,  following  (2),  the  sequence  of  operations  is: 
Choose  an  initial  moment  distribution,  ,  and  put  i  =  1  ; 

2  1  - 

Find  the  curvature  distribution,  c'.'  *  ,  and  put  B._.  =/  ^.t'.'dn  ; 

l  l  i~i  q  li 


Put  t.  -  Bji.ij.  q  -  ; 


1 


Integrate  tV  twice  to  obtain  and  put  a.  =  J  p£.dn  ; 

l  i  i  0  i 

1  1 

Find  the  bending  moment  distribution,  =  /  dHj  J  yt-dri  ; 

1  n  ni  1 

Make  it  orthogonal  to  distributions  in  previous  modes, 

*i+i  =  *i+i  '  Vi  ‘  Bi-i*i-i  *  where  Vo  =  0  5 


>  (9) 


(7)  Add  1  to  i  ,  repeat  from  (2)  until  required  number  of  £.  have  been 
obtained  when  stop  at  (4) . 


All  the  operations  are  numerical,  of  course,  and  the  variables  $,  ij>,  C,  etc 
are  represented  by  their  values  at  a  large  number  of  spanwise  stations.  Again 


the  a.  and  |3.  form  the  leading  diagonal  and  sqb-  and  super-diagonal  of  the 
matrix  A  in  equation  (4).  Similarly  the  eigenvectors  are  given  by 

?a  -  Zqb  (10) 

where  Z  =  (?  ,  $2,  ...,  t^)  ,  and  here  the  are  column  vectors  of  the  values 

of  ?  . 

At  stage  (6)  in  each  sequence  the  orthogonality  to  previous  modes  can 

3 

deteriorate  catastrophically  due  to  cancellation  and  it  is  necessary  that 
j  (in  the  present  case)  be  re-orthogonalised.  This  can  be  done  by  the 
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the  sequence  being  repeated  with  the  same  value  of  j  until 


is  smaller  than 


some  number  e  whose  size  depends  upon  the  numerical  precision  of  the  computer. 
+  )  must  be  reorthogonalised  to  all  Oj  ,  j  =  1,  i  . 

3  SOME  DETAILS  OF  THE  PROGRAM 


The  part-span  integrals  are  evaluated  by  the  trapezium  rule  and  full-span 
integrals  by  Simpson’s  rule. 

The  initial  moment  distribution,  ^  ,  taken  is  that  due  to  the  deflection 
of  the  wing  under  a  combination  of  discrete  heave  load  and  torsional  moment  at 
l  span.  This  is  in  the  spirit  of  the  recommendation  of  Lanczos  (Ref  2,  p.270) 
that  the  initial  vector  should  be  an  arbitrary  vector  multiplied  twice  by  the 
dynamical  matrix. 


The  value  of  the  accuracy  criterion,  e  ,  is  obtained  by  factoring  a  user- 
supplied  value  by  an  approximation  to  the  next  value  of  6  so  that  the  criterion 
is  one  of  relative  accuracy  rather  than  absolute  accuracy. 

Normal  modes  are  calculated  by  using  the  procedure  of  Ref  6. 

Discrete-load  modes  are  included  after  the  normal  modes  of  the  bare  wing 
have  been  calculated.  In  their  initial  form  they  are  the  modes  for  discrete 
force  in  heave  and  moments  in  roll  and  pitch  at  the  loaded  station.  Thus  three 
extra  modes  are  included  for  each  loaded  station.  The  modes  are  subsequently 
made  orthogonal,  with  respect  to  stiffness,  to  each  other  and  to  the  normal 
modes  of  the  bare  wing.  The  orthogonalisation  used  is  that  of  sequence  (11). 

It  is  repeated  because  the  initial  discrete-load  modes  will  be  similar  in  shape 
to  the  bare-wing  modes  already  included.  The  modes  are  normalised  so  that  their 
stiffness  coefficients  are  unity.  The  inertia  matrix  is  transformed  to  tri¬ 
diagonal  by  the  procedure  tred.2  of  Ref  6. 


4  RESULTS  OF  TEST  CALCULATIONS 

The  accuracy  of  an  approximate  normal  mode  is  usually  judged  by  the 
difference  between  its  frequency  and  the  exact  frequency.  When  the  data  for  a 
normal-mode  calculation  are  in  the  form  of  assumed  modes  whose  inertia  and 
stiffness  coefficients  are  known  exactly  the  approximate  frequencies  will  be 
equal  to  or  higher  than  the  exact  frequencies  and  generally  the  percentage  error 
will  increase  motonically  with  the  order  of  the  mode.  This  was  often  the  case  in 
the  present  calculations  but  not  always.  Sometimes  the  error  decreased  after 
first  increasing  and  in  these  cases  the  mode  at  which  the  decrease  first  occurred 
and  higher-order  modes  were  considered  inaccurate.  Another  indication  of 
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numerical  inaccuracy  in  the  formation  of  the  Lanczos  modes  was  the  appearance  of 
spurious  modes  so  that  two  modes  appeared  with  the  same  characteristics  at 
different  frequencies.  Higher-order  modes  (and,  of  course,  the  spurious  modes) 
were  considered  inaccurate. 

4 . I  Bare  wing 

The  majority  of  test  calculations  were  made  for  a  uniform  beam  because  the 
exact  natural  frequencies  are  known.  There  were  two  program  variables  to 
investigate:  one  was  the  number  of  spanwise  elements  used  in  the  numerical 
integration  and  the  other  was  the  number  of  Lanczos  modes  calculated.  The  amount 
of  computer  store  needed  is  roughly  proportional  to  the  product  of  these  numbers 
and  it  is  important  to  know  which  to  reduce  if  storage  is  limited. 

Calculations  were  made  for  a  beam  in  flexure  and  the  results  are  given 
in  Table  I.  Within  the  range  covered  the  accuracy  of  the  modal  frequencies 
increases  as  the  number  of  spanwise  elements  is  increased  but  there  is  a  limit 


Table  1 

Uniform  beam  in  flexure 


No.  of 
integration 
elements 

No.  of  Lanczos 
modes  calculated 

No.  of  normal  modes  found  with 
frequency  error  less  than: 

1% 

0.  1% 

500 

9 

8 

7 

400 

1  1 

9 

8 

300 

15 

12 

7 

9 

8 

7 

200 

22 

6 

5 

15 

10 

5 

100 

22 

7 

2 

15 

8 

2 

9 

7 

2 

after  which  further  increases  do  not  seem  worthwhile.  When  the  number  of 
Lanczos  modes  calculated  is  increased,  the  number  of  modes  found  to  the  higher 
accuracy  does  not  increase.  However,  the  number  of  modes  found  to  the  lower 
accuracy  does  vary  with  the  number  of  Lanczos  modes.  For  300  elements  the 
number  found  to  the  lower  accuracy  increases  with  the  number  of  Lanczos  modes 
but  with  lesser  numbers  of  elements  the  number  of  normal  modes  found  decreases 
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as  the  number  of  Lanczos  modes  is  increased.  In  both  the  200-  and  100-element 
cases,  when  the  number  of  Lanczos  modes  was  22  the  sequence  of  good  inodes  was 
terminated  by  a  spurious  mode.  In  both  cases,  when  the  number  of  Lanczos  modes 
was  15,  the  limit  was  taken  as  the  highest  mode  whose  frequency  error  was 
greater  than  that  of  the  previous  mode.  From  these  results  it  is  more  important 
to  take  a  large  number  of  elements.  In  fact  it  is  dangerous  to  extract  a  large 
number  of  Lanczos  modes  because  of  the  effect  of  cancellation  which  the  results 
for  22  modes  illustrate.  The  number  of  accurate  modes  found  will  obviously 
depend  on  the  numerical  precision  in  the  computer.  All  the  present  calculations 
were  made  using  a  computer  in  which  the  mantissa  is  37  bits  in  length.  A  greater 
number  of  accurate  modes  will  be  found  in  toto  if  there  is  torsion  as  well  as 
flexure.  A  calculation  with  both  flexure  and  torsion,  on  a  beam  described  below, 
gave  15  satisfactory  normal  modes  from  200  elements  and  22  Lanczos  modes. 

Calculations  of  the  torsion  modes  with  the  beam  having  300  elements  gave 
the  results  of  Table  2.  The  modal  displacements  were  within  about  1%  of  the 
exact  values.  To  test  the  program  on  a  non-uniform  beam  the  graver  flexural 

Table  2 

Uniform  beam  in  torsion,  300  elements,  frequency  error  <0.1% 


No.  of  Lanczos  modes 

15 

9 

5 

No.  of  normal  modes 

10 

5 

2 

normal  modes  of  truncated  cones  with  tip  diameters  0.4  and  0.5  of  the  root 
diameter  were  calculated  and  compared  with  the  reference  values  provided  by 
Downs^.  These  calculations  were  made  using  500  elements  and  nine  Lanczos 
modes.  In  each  case  the  first  seven  modes  were  within  0.1%  on  frequency. 

The  fundamental  frequencies  found  were  below  those  given  by  Downs  but  within 

0.01%. 

To  test  the  ability  of  the  program  to  deal  with  modes  close  in  frequency 
the  normal  modes  of  a  uniform  beam  with  both  flexure  and  torsion  allowed  were 
calculated.  There  was  no  coupling  of  flexure  and  torsion  and  the  exact 
frequency  of  the  fundamental  torsion  mode  was  within  0.01%  of  the  exact 
frequency  of  the  first  overtone  flexure  mode.  The  results  of  the  calculation  are 
given  in  Table  3.  The  criterion  for  acceptability  in  this  case  was  the  quality 
of  the  orthogonality  of  flexure  and  torsion.  Because  there  was  no  coupling 
between  flexure  and  torsion  the  normal  modes  calculated  should  have  been  either 
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Table  3 


Uniform  beam  with  close  frequencies  in  flexure  and  torsion, 

300  elements 


No.  of  Lanczos  modes 

15 

13 

I  1 

9 

8 

7 

6 

5 

No.  of  normal  modes 

9 

8 

8 

5 

5 

5 

1 

1 

pure  flexure  or  pure  torsion.  However  calculated  flexural  normal  modes  con¬ 
tained  some  torsion  and  vice  versa  but  the  impurity  could  be  satisfyingly  small. 
There  was  always  a  distinct  division  between  modes  with  satisfactory  and 
unsatisfactory  orthogonality  and  the  frequency  error  of  the  modes  that  were 
acceptable  was  less  than  0.05%.  It  is  noticeable  that  the  number  of  approximate 
normal  modes  does  not  increase  regularly  as  the  number  of  Lanczos  modes  is 
increased.  In  this  example  not  only  are  the  second  and  third  frequencies  almost 
the  same  but  the  fourth  and  fifth  frequencies  are  within  7%  of  each  other  and  the 
sixth  and  seventh  frequencies  are  within  10%  of  each  other  and  so  the  observation 
of  Lanczos  on  systems  with  widely-spaced  frequencies  to  the  effect  that  another 
normal  mode  is  obtained  for  every  intermediate  mode  added  to  the  calculation  does 
not  apply. 

4.2  Loaded  wing 

Test  calculations  were  made  for  a  loaded  wing  to  prove  that  part  of  the 
program  in  which  discrete-load  modes  are  involved .  In  the  tests  very  large 
masses  and  moments  of  inertia  were  put  at  the  tip  so  that  the  overtone  modes 
could  be  compared  with  the  modes  of  a  propped  cantilever  and  built-in  beam. 

In  all  the  calculations  200  elements,  eight  Lanczos  modes  and  five  bare-wing 
normal  modes  were  used. 

For  only  a  mass  at  the  tip  (propped  cantilever)  the  first  five  frequencies 
of  the  loaded  wing  were  found  to  within  0.1%.  For  both  mass  and  rolling  inertia 
at  the  tip  (built-in  beam)  the  first  four  frequencies  were  found  to  within  1%. 

For  only  pitching  inertia  at  the  tip  (built-in  beam)  the  first  four  torsional 
frequencies  were  found  to  within  0.1%.  These  test  cases  are  severe  and  show 
that  good  results  can  be  obtained  even  when  the  numbers  of  integration  elements 
and  bare-wing  modes  used  are  not  large. 

4.3  Comments  on  use 

From  the  results  obtained  the  method  should  give  good  approximations  to  as 
many  of  the  graver  normal  modes  as  are  needed  in  aeroelastic  investigations. 
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In  view  of  the  results  in  Table  3  it  is  probably  wise  to  use  at  least  seven  inter¬ 
mediate  modes  but  it  is  quite  possible  that  five  good  normal  modes  can  be  obtained 
from  these  seven  intermediate  modes  if  the  beam  is  divided  up  into  an  adequate 
number  of  elements.  If  five  or  more  modes  are  needed  however  it  is  better  to  use 
more  intermediate  modes  but  there  seems  to  be  little  point  in  using  even  as  many 
as  twice  the  number  of  normal  modes  required.  Indeed,  if  the  number  of  normal 
modes  required  is  large,  some  restraint  needs  to  be  shown  in  the  number  of  inter¬ 
mediate  modes  used  if  numerical  difficulties  are  to  be  avoided  but  the  position 
is  eased  if  the  beam  can  distort  in  both  flexure  and  torsion. 

Little  mention  has  been  made  of  the  accuracy  of  the  shape  of  the  calculated 
mode  because  evaluating  it  is  difficult.  However,  the  calculated  modes  were 
studied  and  it  is  thought  that  the  maximum  error  in  deflection  will  be  of  the 
order  of  1%  for  a  mode  with  a  frequency  error  of  less  than  0.1%. 

5  CONCLUSIONS 

The  Lanczos  method  of  minimised  iterations  has  been  developed  to  calculate 
the  graver  normal  modes  of  cantilever  beam-like  structures  directly  from  the 
distributions  of  their  mass  and  stiffness.  The  method  is  not  one  for  purely 
mechanical  application  but  the  computer  program  written  has  given  accurate 
results  for  as  many  modes  as  will  be  required  in  most  practical  cases.  No 
definite  rules  are  given  for  the  number  of  intermediate  modes  that  should  be 
used  because  this  will  depend  on  the  numerical  precision  of  the  computer  being 
used.  However  illustrative  results  have  been  found  which  indicate  the  relative 
merits  of  having  a  large  number  of  integration  elements  and  a  large  number  of 
intermediate  modes.  The  program  also  includes  provision  for  discrete-load 
modes  and  the  application  of  these  has  been  demonstrated  by  calculating  the 
normal  modes  of  beams  carrying  large  concentrated  masses. 
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